Lecture 14: Spatial GLMs

Class Intro

Intro Questions

  • Describe geometric anisotropy. Discuss this in the context of a data set.

  • For Today:
    • Finish geometric anisotropy
    • Spatial Generalized Linear Models

Geometric Anisotropy

Geometric Anisotropy Model

  • Let \(\boldsymbol{Y}(\boldsymbol{s}) = \mu(\boldsymbol{s}) + w(\boldsymbol{s}) + \epsilon(\boldsymbol{s})\), thus \(\boldsymbol{Y}(\boldsymbol{s}) \sim N(\mu(\boldsymbol{s}), \Sigma(\tau^2, \sigma^2, \phi, B))\), where \(B = L^T L\).
  • The covariance matrix is defined as \(\Sigma(\tau^2, \sigma^2, \phi, B)) = \tau^2 I + \sigma^2 H((\boldsymbol{h}^T B \boldsymbol{h}^T)^{\frac{1}{2}}),\) where \(H((\boldsymbol{h}^T B \boldsymbol{h}^T)^{\frac{1}{2}})\) has entries of \(\rho((\boldsymbol{h_{ij}}^T B \boldsymbol{h_{ij}}^T)^{\frac{1}{2}}))\) with \(\rho()\) being a valid correation function, typically including \(\phi\) and \(\boldsymbol{h_{ij}} = \boldsymbol{s_i} - \boldsymbol{s_j}\).

Geometric Anisotropy Visual

  • Consider four points positioned on a unit circle.

  • How far apart are each set of points?

Geometric Anisotropy Exercise 1

Now consider a set of correlation functions. For each, calculate the correlation matrix and discuss the impact of \(B\) on the correlation. Furthermore, how does B change the geometry of the correlation?

  1. \(\rho() = \exp(-\boldsymbol{h_{ij}}^T B \boldsymbol{h_{ij}}^T)^{\frac{1}{2}})),\) where \(B = \begin{pmatrix} 1 & 0 \\ 0 & 1 \\ \end{pmatrix}\)

  2. \(\rho() = \exp(-\boldsymbol{h_{ij}}^T B \boldsymbol{h_{ij}}^T)^{\frac{1}{2}})),\) where \(B = \begin{pmatrix} 2 & 0 \\ 0 & 1 \\ \end{pmatrix}\)

  3. \(\rho() = \exp(-\boldsymbol{h_{ij}}^T B \boldsymbol{h_{ij}}^T)^{\frac{1}{2}})),\) where \(B = \begin{pmatrix} 3 & 1 \\ 1 & 1 \\ \end{pmatrix}\)

Geometric Anisotropy: Solution 1

\(\rho() = \exp(-\boldsymbol{h_{ij}}^T I \boldsymbol{h_{ij}}^T)^{\frac{1}{2}}))\)

Implied Distance

0.00 1.41 1.41 2.00
1.41 0.00 2.00 1.41
1.41 2.00 0.00 1.41
2.00 1.41 1.41 0.00

Correlation

1.000 0.243 0.243 0.135
0.243 1.000 0.135 0.243
0.243 0.135 1.000 0.243
0.135 0.243 0.243 1.000

Geometric Anisotropy: Solution 2

\(\rho() = \exp(-\boldsymbol{h_{ij}}^T B \boldsymbol{h_{ij}}^T)^{\frac{1}{2}})),\) where \(B = \begin{pmatrix} 2 & 0 \\ 0 & 1 \\ \end{pmatrix}\)

Implied Distance

0.00 1.73 1.73 2.83
1.73 0.00 2.00 1.73
1.73 2.00 0.00 1.73
2.83 1.73 1.73 0.00

Correlation

1.000 0.177 0.177 0.059
0.177 1.000 0.135 0.177
0.177 0.135 1.000 0.177
0.059 0.177 0.177 1.000

Geometric Anisotrop: Solution 3

\(\rho() = \exp(-\boldsymbol{h_{ij}}^T B \boldsymbol{h_{ij}}^T)^{\frac{1}{2}})),\) where \(B = \begin{pmatrix} 3 & 1 \\ 1 & 1 \\ \end{pmatrix}\)

Implied Distance

0.00 1.41 2.45 3.46
1.41 0.00 2.00 2.45
2.45 2.00 0.00 1.41
3.46 2.45 1.41 0.00

Correlation

1.000 0.243 0.086 0.031
0.243 1.000 0.135 0.086
0.086 0.135 1.000 0.243
0.031 0.086 0.243 1.000

More Geometric Anisotropy

  • The matrix \(B\) relates to the orientation of a transformed ellipse.
  • The (effective) range for any angle \(\eta\) is determined by the equation \[\rho(r_\eta(\tilde{\boldsymbol{h}}_{\eta}^T B \tilde{\boldsymbol{h}}_{\eta}^T)^{\frac{1}{2}}) = .05,\] where \(\tilde{\boldsymbol{h}}_{\eta}\) is a unit vector in the direction \(\eta\).

Fitting Geometric Anisotropy Models

  • Okay, so if we suspect that geometric anisotrophy is present, how do we fit the model? That is, what is necessary in estimating this model?
  • In addition to \(\sigma^2\) and \(\tau^2\) we need to fit \(B\).
  • What about \(\phi\)? What is \(\phi\) when \(B = \begin{pmatrix} 1 & 0 \\ 0 & 1 \\ \end{pmatrix}\)?

Priors for B

  • While \(B\) is a matrix, it is just another unknown parameter.
  • Hence, to fit a Bayesian model we need a prior distribution for \(B\).
  • One option for the positive definite matrix is the Wishart distribution, which is a bit like a matrix-variate gamma distribution.

Spatial GLMS

Data Viz Motivation

Ozone Exceedance in Colorado

Generalized Linear Model Notation

There are three components to a generalized linear model:

  1. Sampling Distribution: such as Poisson or Binomial
  2. Linear combination of predictors: \(\eta = X\beta\)
  3. A link function to map the linear combination of predictors to the support of the sampling distribution.

Logistic Regression Overview

Write out the complete model specification for logistic regression.

  • Assume \(Y_i\) is the binary response for the \(i^{th}\) observation, \[\begin{eqnarray*} Y_i &\sim& Bernoulli(\pi_i)\\ logit(\pi_i) &=& X_i \beta, \end{eqnarray*}\]
  • where \(logit(\pi_i) = log \left(\frac{\pi_i}{1-\pi_i}\right)\)

glm()

Interpret this output

CO <- CO %>% mutate(north = as.numeric(Latitude > 38 ))
glm(Exceedance~north, family=binomial(link = 'logit'),data=CO) %>% summary
## 
## Call:
## glm(formula = Exceedance ~ north, family = binomial(link = "logit"), 
##     data = CO)
## 
## Deviance Residuals: 
##    Min      1Q  Median      3Q     Max  
## -2.313   0.378   0.378   0.378   1.177  
## 
## Coefficients:
##               Estimate Std. Error z value Pr(>|z|)  
## (Intercept) -1.695e-15  8.165e-01   0.000   1.0000  
## north        2.603e+00  1.097e+00   2.372   0.0177 *
## ---
## Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
## 
## (Dispersion parameter for binomial family taken to be 1)
## 
##     Null deviance: 28.708  on 34  degrees of freedom
## Residual deviance: 22.873  on 33  degrees of freedom
## AIC: 26.873
## 
## Number of Fisher Scoring iterations: 5

Logistic Regression in JAGS

JAGS specification for logistic regression

model{
  # likelihood
    for (i in 1:N){
        y[i] ~ dbern(p[i])
        p[i] <- 1 / (1 + exp(-eta[i]))
        eta[i] <- int + x[i] * north
    }
    # prior
    int ~ dnorm(0, 1E-12)
    north ~ dnorm(0, 1E-12)
}

Spatial Logistic Regression

  • Assume \(Y(\boldsymbol{s_i})\) is the binary response for \(\boldsymbol{s_i}\), \[\begin{eqnarray*} Y(\boldsymbol{s_i})|\beta, w(\boldsymbol{s_i}) &\sim& Bernoulli(\pi(\boldsymbol{s_i}))\\ logit(\pi(\boldsymbol{s_i})) &=& X(\boldsymbol{s_i}) \beta + w(\boldsymbol{s_i}), \\ \end{eqnarray*}\]
  • where \(\boldsymbol{W} \sim N(\boldsymbol{0},\sigma^2 H(\phi))\)

Spatial Logistic Regression in JAGS

logistic.spatial <- "model {
  for (i in 1:N) {
      Y[i] ~ dbern(p[i])
      logit(p[i]) <- mu[i] + w[i]
      mu[i] <- beta[1]+beta[2]*x.north[i]
    muW[i] <- 0
  }
  # process
  w[1:N] ~ dmnorm(muW[],Omega[,])

  # priors
  beta[1] ~ dnorm(0.0,1E-8)
  beta[2] ~ dnorm(0.0,1E-8)
  spat.prec ~ dgamma(0.001, 0.001)
  sigmasq <- 1/spat.prec
  phi ~ dunif(1.5,2.5)  
    
  #build omega
  for (i in 1:N){
    for (j in 1:N){
      H[i,j] <- (1/ spat.prec) * exp(-phi * d[i,j])
    }
  }
  Omega[1:N,1:N] <- inverse(H[1:N,1:N])
}"

spGLM()

m.1 <- spGLM(y~1, 
             family="binomial", 
             coords=coords, 
             weights=weights, 
             starting=list("beta"=beta.starting, "phi"=0.06,"sigma.sq"=1, "w"=0),
             tuning=list("beta"=beta.tuning, "phi"=0.5, "sigma.sq"=0.5, "w"=0.5),
             priors=list("beta.Normal"=list(0,10), "phi.Unif"=c(0.03, 0.3), "sigma.sq.IG"=c(2, 1)),
             amcmc=list("n.batch"=n.batch, "batch.length"=batch.length, "accept.rate"=0.43),
             cov.model="exponential", 
             verbose=TRUE, 
             n.report=10)

Spatial Poisson Regression

Motivation

Poisson Regression Overview

Write out the complete model specification for Poisson regression.

  • Assume \(Y_i\) is the count response for the \(i^{th}\) observation, \[\begin{eqnarray*} Y_i &\sim& Poisson(\lambda_i)\\ \log(\lambda_i) &=& X_i \beta, \end{eqnarray*}\]
  • thus \(\exp(X_i \beta) \geq 0\)
  • next write out a Poisson regression model with spatial random effects

Spatial Poisson Regression

  • Assume \(Y_i\) is the count response for \(\boldsymbol{s_i}\), \[\begin{eqnarray*} Y(\boldsymbol{s_i}) &\sim& Poisson(\lambda(\boldsymbol{s_i}))\\ \log(\lambda(\boldsymbol{s_i})) &=& X(\boldsymbol{s_i}) \beta + w(\boldsymbol{s_i}), \end{eqnarray*}\]
  • where \(\boldsymbol{W} \sim N(\boldsymbol{0},\sigma^2 H(\phi))\)

Fitting Spatial Poisson Regression Models

  1. spGLM() accepts binomial and Poisson data,
  2. JAGS can also be used to fit spatial GLMS models, or
  3. Try writing a sampler from scratch.